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Abstract 

The aim of this study is to find a generic method for generating a path 
of the solution of a given stochastic differential equation which is more 
efficient than the standard Euler-Maruyama scheme with Gaussian incre¬ 
ments. First we characterize the asymptotic distribution of pathwise error 
in the Euler-Maruyama scheme with a general partition of time interval 
and then, show that the error is reduced by a factor {d + 2)/d when using 
a partition associated with the hitting times of sphere for the driving d- 
dimensional Brownian motion. This reduction ratio is the best possible in 
a symmetric class of partitions. Next we show that a reduction which is 
close to the best possible is achieved by using the hitting time of a moving 
sphere which is easier to implement. 


1 Introduction 

Various stochastic phenomena have been modeled in terms of the solution 
X = (X^,..., XP) of a stochastic differential equation (SDE) 

d 

dX'(s) = E fj{X{s))dWi{s) (1.1) 

i=Q 

on a domain D c IRP, where W = (W^,..., W‘^) is a d-dimensional standard 
Brownian motion, W°(s) = s is the time coordinate, and / = {/p 1 < i < p,0 < j < 

d} : D ^ IRP ® is a continuously differentiable fimction. The Monte Carlo 

'Department of Mathematics, Osaka University, 1-1 Machikaneyama, Toyonaka, Osaka, Japan, 
email: fukasawaOmath. sci .Osaka-u.ac. jp, web: www.math. sci .osaka-u. ac.jp/-fukasawa 
'''Mathematical Institute, University of Oxford, ROQ, Woodstock Road, Oxford 0X2 6GG, UK. 
email: jan.obloj@maths.ox.ac.uk, web: www.maths.ox.ac.uk/people/jan.obloj Jan Obtoj 
gratefully acknowledges funding received from the European Research Council under the Euro¬ 
pean Union's Seventh Framework Programme {EP7/2007-2013) / ERG grant agreement no. 335421 
and is also thankful to the Oxford-Man Institute of Quantitative Finance and St John's College in 
Oxford for their financial support. 


1 



simulation is a powerful and very popular approach to study such a stochastic 
model. The standard method for generating a path which follows the SDE 
is the Euler-Maruyama scheme that constructs an approximating sequence of 
processes X” = (X”'\..., X”'P) as 

d 

X"(0) = X(0), X"''(S) = X”'‘(7i”) + fiQC\nl)){W{s) - W{ul)) (1.2) 

j=0 

for s e (tt",, tt” _^j], where tt" = {nj} is an increasing sequence of stopping times, 
usually chosen to be a deterministic sequence such as tt” = m/n. The variable n 
controls the computational effort of this construction. We naturally expect that 
X” —> X in some sense as M —> oo. The attractive features of the Euler-Maruyama 
scheme include its validity under degenerate diffusion coefficients with mild 
regularity, intuitive construction and easy implementation. See Kloeden and 
Platen [TT] for some elementary properties of this and other related methods. 
Kurtz and Protter ca studied the limit of the approximation error process 

W = {W'\U"'‘{s) = Vw(X”'*(s) - X‘(s)) (1.3) 

as n —> oo. They showed that if the sequence of {d + l)^-dimensional processes 
Z” = 0<l<d,0<i<d] defined by 

Z"'''^(f) = (W'(s) - W'«))dW^(s) (1.4) 

m=0 

is "good" and converging to a semimartmgale Z = {Z^'l} in law, then Lf" con¬ 
verges in law toU = (W, LfP), the solution of 

dU\s) = Y dkfj{X{s))U^{s)dWi{s) - Y dkfj{X{s))fl^{X{s))dZ‘'i{s). (1.5) 
i,k i,kJ 

Since this SDE for U is affine, we may write LI in a more explicit form. In 
particular if (Z^'l, W') = 0 for all I, j, i, then 

U{t) = Y{t) f Y(s)-^dF(s), dF\s) ^ -Y^kfl(X(s))fiHX(s))dZ%), (1.6) 

;,kj 

where Y = {Y”'*’; 1 < a, b < p} is the solution of 

dY"''’(s) = Y <?ic//(X(s))Y^'^(s)dW^'(s), Y”'''(0) = (1.7) 

j.k 


where 6“'^ = is the Kronecker's delta. For example in the case of p = 1, 
U{t) = -Y{t)Y f Y(s)-^f'(X(smX(s))dZ%), 

I tn do 
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and 



It is important to note that Y does not depend on tt” . In consequence, the limit 
distribution U of U" is determined by the limit distribution Z of Z" in a linear 
maimer. In the equidistant case n” = min, as shown by Kurz and Protter flZl 
and Jacod and Protter ||9l, we have = Zi’^ = 0 for all 0 < / < d and 

^f2Z‘'i, l<l,i<d 

is a d^-dimensional Brownian motion independent of W. In parficular, we 
have (Z^'i, W') = 0 for all I, j, i. and that the distribution of U{t) is condifionally 
Gaussian. 

In this paper we consider more general sequences of partifions { tt" }. Our first 
main result states a central limit theorem for the error process and provides a 
convenient characterization of the resulting process Z. This allows us to study 
the efficiency of different sequences of partitions. Our second main results 
establishes a uniform lower bound on the expected asymptotic error U and 
shows that the bound is attained by partitions associated with the hitting times 
of a sphere for W. Note that in such a scheme both the time steps Um - Hm-i and 
the W increments, j are random. The later are uniformly distributed 

on a sphere and in particular in the one-dimensional case take just two values. 
In contrast, in the equidistant case the time steps are deterministic (and equal) 
and all the randomness is coming from the increments of W. Newton IT^ 
and Fukasawa m studied the hitting time scheme in the one-dimensional case. 
Cambanis and Hu HI and Miiller-Gronbach m gave optimality results among 
irregular and adaptive schemes which still use Gaussian increments. Our 
framework of discretisation is more general than these preceding studies. Our 
result is closely related to a recent work by Landon flSlI , where a sequence 
of hitting times of ellipsoids is derived as an asymptotically optimal scheme. 
In this paper, we restrict schemes to be symmetric in a certain sense because, 
among other reasons, there is unlikely to be a realistic computational algorithm 
to implement asymmetric schemes. Our framework therefore excludes hifting 
times of ellipsoids (except spheres). 

Our theoretical analysis of efficiency described above does not take into ac¬ 
count the complexity of simulating random variables with a given distribution, 
which may be challenging for hitting times of spheres in higher dimensions. 
We argue that in practice one should look for a scheme where both the time 
steps and the spatial increments are random and all are easily generated to¬ 
gether. We achieve this adapting the moving sphere approach in Deaconu and 
Herrmann 0. This scheme is easier to implement and enjoys a better accuracy 
than the standard Gaussian scheme. In fact, it modifies the standard one only 
by replacing 
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with 


^m-1 ” itlt ■/ 

U^ml 

where 

,7 _ + 2Em 

^m ~ j 

a 

and Em ~ Exp(l) and Nm ~ NiO, Id) are independent iid sequences. It improves 
the accuracy of the Monte Carlo simulation even after taking into account a 
slight increase of computational time due to the one additional generation of 
random variable and the calculation of the exponential function each step. 
It also has the further advantage that both the time and W increments are 
bounded. This enables us to control the size of each increments of X" and deal 
with SDE on a bounded domain, see Milstein and Tretyakov IIT4l . or devise 
efficient pricing of path dependent options e.g. barrier options. 

This paper is organised as follows. In Section 2 we present a central limit 
theorem for discretisation error. In Section 3 we study several examples of 
schemes and discuss their effectiveness. In particular, we show one of them to 
be attractive in terms of both error magnitude and computational costs. 


2 Central limit theorem for the error process 

Here we present a central limit theorem for the asymptotic error process !T" = 
Vn(X" - X). Let (Q, P) be a probability space. We suppose that the SDE 
O admits a unique strong solution X which does not explode and remains 
in a given open cormected domain D c PP. We further assume there exists 
a sequence of compact sets with each K,„ being a subset of the interior of 
K^+i such that = D and “ as m —> 00 , where 

Tm = inf(s > 0;Xf ^ Km}- 


Denote by (;?s}s>o the augmentation of the natural filtration generated by W. A 
partition n = {7im}m>o is a sequence of increasing stopping times with tiq = 0 
and limm^oo nm = Eor a partition n, put 


— Tim TIm-1, 


A” W = W„ - W.„ 

T"" = 'F 

m ■' 

N" = min{m > 0; Um > z}. 


Note that AT" is the number of discretisation steps required to generate a path 
up to a finite stopping time z. In this section we do not take the computational 
difficulty to generate Amii and A^W into accoimt. Hence we take N" as a 
measure of computational effort associated with the partition n. Notice that N” 
is a stopping time with respect to the discrete filtration 
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For a given sequence of partitions n", we denote by A” W, N" and 
the corresponding quantities A^^W, N"™ and for brevity. We also let 
7t"(s) = 7t” for s e [nlJ,, For example if {n"} is the equidistant scheme, 

that is, TiJJ, = m/n, then we have N" = [nx] ~ mt and 7t”(s) = [nsj/n ~ s as 
n ^ c)o. To see what happens if {n”} is stochastic, let us first consider a simple 
example which is known as the adaptive scheme, see IflSl . Let G be a positive 
continuous adapted process and define ttJJ, by 


< = 0, 7l” = 7 t” + 


nG{nl) 


( 2 . 1 ) 


Then 


N" ^ 

-A = G(7t” _i)A„,h” ~ I 

” ^1 -'o 


G(s)ds 


as n ^ oo. 

More generally if there exist adapted processes G" and G such that 


ie[a,„+i7t”|;t;”] 


1 

nG"«) 


( 2 . 2 ) 


and 




sup |G”(7 i"(s)) - G(s)| 0, 2^ E [|A„,7i”niF;”_J ^ 0 (2.3) 


0<S<T 


m=l 


in probability, then we have 


r 

n Jo 


G(s)ds 


(2.4) 


in probability as n —> oo by a simple application of fhe Lenglarf inequality 
(see e.g. Lemma A.2 in Fukasawa [Zl). The computational effort is therefore 
controlled by the process G. We want to consider here schemes {n’J which 
satisfy (12.211 and (12.3b together with a mild symmetry requirement: 




where L' = 


f 


(2.5) 


(W\s) - W\nl_,))iWi(s) - Wi(nl_,))ds, 


for each n, m and 1 < i, j < d with i + j, and further 


sup |H"( 7 t"(s)) - H(s)| -> 0, 

0<S<T 

N? n; (2.6) 

2^ E[|A,7t"n!r”_J ^0, 0 

m=l m=l 
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in probability as n —> oo, where H" and H are some adapted processes. Note 
that the above conditions hold with G„ = G and H„ = H = 3/G if n" is given by 
with a continuous adapted process G. We can treat a scheme like 


71" =0, 


n(G«) A n) 


as well by taking G" = G A n. 

We now state our central limit theorem for the error process using discretisa¬ 
tion schemes as above. Our standing assumption is that the Euler-Maruyama 
approximation X" is well-defined by II1.2II , i.e., X"(7t|5,) keeps inside the domain 
of / up to certain time t > 0. Fix such t > 0. The error process Li" is then 
well-defined by l ll.3ll up to the time t. 


Theorem 2.1 Let {rj^} be an increasing sequence of stopping times with 

lim P[?]a: < t] = 0. 

k—*oo 

Let {tt"} be a sequence of partitions satisff/ing,for all k, (12.2P , (12.3P , ( 12.5P and (12.6P with 
locally integrable processes G and Li and t = ^f supg<<,<j |!J"(s)| is tight, then the 

C{[0,t],W)-valued random sequence U" converges in lawlo the solution U of il.Sh 
where Z is given by Z^'i = Zi'^ = Z^’^ = 0, 

Z'4)=^ r V^dW"^(s), (2.7) 

V6 Jo 

and W = {W'd is a d^-dimensional standard Brownian motion independent ofW. In 
particidar (I2.6P holds. 


The proof of this result is given in Section 4. We note that the symmetry 
condition (12.5P is not essential and the above result could be extended to asym¬ 
metric schemes. However, we chose to concentrate on the given framework 
for number of reasons. First, asymmetric schemes are in general harder to 
implement. Second, whether or not an asymmetric scheme is superior to the 
standard equidistant one depends on the coefficients of the SDE (|l.lb which 
goes against our main purpose of finding a generic method which uniformly 
improves on the standard method. Third, an asymmetric scheme induces an 
asymptotic bias and an additional source of randomness in the limit, which 
is not preferable from practical point of view. This bias can be corrected but 
the correction makes implementation more complicated. See Fukasawa [3 for 
related results in the one-dimensional case. 


3 Asymptotically efficient schemes 

We apply now Theorem 12.11 to compare and study different discretisation 
schemes. We establish a lower bound on the expected squared asymptotic 
error E[|Lrfp] and exhibit an efficient scheme which attains the bound. We 
start however with the usual benchmark given by the classical schemes with 
Gaussian increments W. 
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3.1 Gaussian schemes 

Consider first {n"} defined by (12.1b with a positive continuous adapted process 
G. As already mentioned, Il2.2b . Il2.3b . (12.5b and (12.6b are satisfied with G„ = G 
and H„ = H = 3/G for any finite stopping time t. By the standard theory of 
the Euler-Maruyama scheme (see e.g. Kloeden and Platen IfTTI l if, say, G = 1, 
D = ]RP and the coefficient / is of linear growth, then 

sup]E[sup |Lt"(s)|^] < oo (3.1) 

n 0<s<t 


for each t > 0, which implies that supg^g^j |Lt"(s)| is tight. By localization 
argument we can easily conclude this tightness without the restriction G = 1. 
Then we can apply Theorem 12.II to have 


U{t) = Y{t) 


f 


y(s)-Mf(s), dF(s) = -— 

V2 y ;t=l 


Z, dkf.{X{s))fHX{s)) 


^/W) 


dW‘'i(s) 


(3.2) 

as the limit of Lt”, where Y = {Y”'*’; 1 < a,b < p] is the solution of (11.7b . 

Next consider n", = g(m/n), where g is a continuously differentiable increas¬ 
ing fimction with g(0) = 0. Since nAn’^ ~ g'{g~^{n’^_^)) as n —> oo, (I2.2b , (12.3b . 
(12.5b and (12.6b are satisfied with G"(s) = l/(g(g“Hs) + 1/w) - s), G = 1/g' o g~^, 
H" = 3/G” and H = 3/G. This scheme has essentially the same structure as (12.1b 
and the tightness of supg<j,<j |G"(s)| is verified by the same manner. The limit of 
Lt" is also given by (13.2b . 

Both of these schemes use conditionally Gaussian increments A” W. They 
are relatively easy to implement and widely used in practice. The particular case 
G = 1 is the standard equidistant scheme. However, when the computational 
effort is measured by (12.4b , they turn out to be inefficient in terms of asymptotic 
error magnitude. 


3.2 Efficient scheme 

By Theorem l2.ll the asymptotic error distribution depends on jn"} only via H 
in (12.6b . while the computational effort in (12.4b does only via G. The smaller 
the H, the smaller the error. More precisely, tor a discretisation scheme which 
satisfies the assumptions in Theorem 12 .1 1 and (13.1b , we have 

f ©(t,s)H(s)ds 
0 

with 0(t, s) > 0 defined purely in terms of {Xs}o<s<f, where we used (11.6b and the 
independence of W and W. In the following theorem we give a lower bound 
on H, for a given G, and exhibit a scheme which attains it. 

Theorem 3.1 Let t > 0 and a locally integrable adapted process G be given. Let {rjk] 
be an increasing sequence of stopping times with limjt^oo Plqt < t] = 0. Let {tt”} he 


E [n|X"(t) - X(t)p] = E [|;J”(t)p] —» E 


7 






a sequence of partitions satisfying (12.2P , ||2.3I I, (I2.5I > and (I2.6I > a locally integrable 
process H and t = pkfor all k. Then, 


Hu> 


3d 


du ® dP a.e. on [0, t] x O. 


(3.3) 


^0-"' 'Wl - ^ ( ’ (3-4) 

satisfies (12.2P , (12.3P , (I2.5P and (l2.6P for t = fjk and attains the equality in (13.3t , where 


{d + 2)Gu 

If G is positive and continuous on [0, t], then the sequence {n”] defined as 
K = 0, n" = inf If > n” : \W{t) - W«)|^ = ^ 


fjk = inf(s > 0 : G(s) > kor G(s) < 1/k}. 

Remark 3.2 We note that the different discretisation schemes we consider mirror 
different pathwise constructions of the stochastic ltd integral. The equidistant scheme 
in (I2lt with G ^ 1 (or other deterministic function) is akin to the approximation of 
stochastic integral discussed in Fdllmer ®. The asymptotically efficient scheme in 
iu in contrast corresponds to discretising the path along "Lebesgue type partition", 
as in Vovk and Perkowski and Prdmel see also Davis et al. ^for a discussion 

and more general partitions. 

Compared with the standard equidistant scheme, the error distribution for the 
scheme ll3.4t uniformly shrinks ; denoting by Ugia and Licauss the limits of Lf" 
associated with (13.4b and (12.lb respectively, we conclude 


Uem'= 


Both ll3.4b and (12.1b require the same computational effort (12.4b . This refines a 
result for one-dimensional case given in Fukasawa 13 . 

To prove the lower bound (13.3b . it is sufficient to establish it for processes 
H” and G”. Given ll2.2b and (12.5b . this is equivalent to 


E [\Al,,Wfy\r:] > (E , Vn, m. 


It follows that Theorem l3.1l is a direct consequence of the following lemma. 
Lemma 3.3 Let a > 0 and 

d 

Q(v) = Y,\W(vt. 

y=i 


Then 


min{E[Q(T)] : t is a stopping time with E[t] = a} = 
and the minimum is attained by 

T = inf |t > 0 : IW(f)|^ = dflj. 


3d^a^ 
d + 2 


(3.5) 

(3.6) 
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The proof of the Lemma is deferred to Section 4. We discuss now how to 
implement (13.41 . Since |W| and W/|W| are independent, conditionally on 
A",W is independent of A„n” and uniformly distributed on the sphere with 
radius 

By a scaling property. 


A” W- 


N 


’ nG(n"_^) INI' 




d 




Tl 


where N ~ N{0, Id) and n is defined by (13.61 with a = 1/d, which has the same 
law as the hitting time of 1 for the d-dimensional Bessel process starting from 
0. Generating A"jW is well discussed and for A„,7 t", it suffices to develop an 
efficient algorithm for generating ti by, say, the acceptance-rejection method. 
An explicit form of the distribution function of ti is given by Ciesielski and 
Taylor ||3). The implementation effort and complexity vary with d. We do not 
pursue this further here. Instead, we provide an attractive alternative in the 
next section. 


3.3 Moving sphere scheme 

We adapt here the moving scheme algorithm presented by Deaconu and Her- 
rmarm 0. The idea is to consider partitions defined by hitting times of a 
sphere with a radius shrinking in time. The rate at which radius shrinks is 
adjusted in such a way that both the time step and the spatial increment have 
explicit distributions which are easy to simulate numerically. Both distribu¬ 
tions are non-trivial in the sense that they admit density on some set of positive 
Lebesgue measure. This is in contrast to the two extreme schemes: the classical 
equidistant scheme in which time steps are deterministic and the asymptotically 
efficient scheme of Theorem 13.11 in which the increment A",W is concentrated 
on a (fixed) sphere. 

Let G be a positive continuous adapted process on [0, f]. We define the 
sequence of partitions {ti"} via 

< = 0, nl,, = inf |s > <;|W(s) - W(h;;)P > G” (3.7) 

where 

1 a ! 

S's) ■ 

Since ip(a) = 0, A„,7t” is bounded by flG|5,. Although ip(Q) = 0, we can show 
A„,7 t" > 0 a.s. by the law of iterated logarithm. Since |W| and W/|W| are 
independent, conditionally on 

(A^n\AlW) ~ 
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where N ~ N{0, Id) and 


= inf{s > 0;|W(s)P > i/»(s)}. 


Now we show that generating a random variable with the same distribution as 
Tij, is quite easy. By Proposition 2 of Deaconu and Herrmarm ||5l, the density of 
is given by 


s i-> 


T{d/2)2'i/V/h 


dt log ■ 


d/2 


Remark that a method of Chen et al. HI can be applied to prove this with a slight 
modification. As shown by Proposition A.l of Deaconu and Herrmarm ||3, we 
have then that 

~ ae~^, 


where Z is a random variable which follows the Gamma distribution with 
shape 1 + d/2 and scale 2/d. Since |Np is independent of N/\N\ and follows the 
Gamma distribution with shape d/2 and scale 2, we can use |Np to generate Z 
as 

Z~ i(|N|2 + 2E), 

where E is an exponentially distributed random variable with mean 1 which is 
independent of N. Thus we have 


A" W) ~ [cUae-^ ^Gl_,adZe-^^^ 


conditionally on . 

Now we show that i2.2L Il2.3b . (12.5b and (12.6b are satisfied. Observe that 


]E[A™^i7i"r”] = aGlEle-^] = Gl = 
]E[|A,H''h!F^_J = a^\GlMe-'^^) = 0{n-% 


-E[Nn = E 
n 


N? 




m=l 



- N" 

<KE 

A^tt" 


m=l 


< K{t + K/n), 


where 

T = rjK = t A inf(s > 0 : G(s) > K or G(s) < 1/fC} 
for any X e N. It follows then 


m=l 
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in probability for all fc > 2 and z = r]x- Further by the same argument as in the 
proof of Lemma l33l we have that 


E[|A” 


/=1 

3|G”|2 


= 


'm+l 




d{d + 2) 'I" ^ d + 2 

3\G-J{d + 2Y^^ _ 1 H«) 

d^l'^{d + fi2 G{nli)' 


d(d + 2) 

fl^E[Z2e-2z 


E[|A" wn;?;" 


where 



r(d) = 


{d + 2)'^+2 

d‘</2(d + 4)(‘^+4)/2 • 


(3.8) 


Figure [T] plots the reduction ratio r(d) with efficiency bound d/(d + 2) from 
Il3.3b in red. As is clearly seen, r{d) < 1, which means (13.7b is superior to 
Gaussian schemes. We can show that r(d) < dl{d + 1), which implies that the 
moving sphere scheme keeps the advantage even after taking into account that 
the proposed method requires generating one additional exponential random 
variable in each step. 

Further from Figure [TJ we find fhat the reduction ratio r{d) is close to the 
best possible value d/{d + 2) attained by ll3.4b . Taking computational costs for 
generating increments A^n” into accoimt, (13.7b is likely to be more efficient than 
ll3.4b for most applications. 

Finally, the moving sphere scheme has an additional advantage that both 
the increments AmTi” and A” W are bounded. More precisely 


A,„7t"<flG” and |W(s) - W«)| <-G” (3.9) 

for all m and n. In consequence, by changing G in an adapfed way, we can 
control the size of each increments of X". This is necessary fo deal with the 
SDE on a bounded domain, see Deaconu and Herrmann or Milstein and 
Tretyakov 1141 . Further, this also allows to adapt the scheme to obtain a greater 
accuracy of certain path information. Consider for example pricing of a barrier 
option - we may monitor the barrier crossing of our simulafed pafhs with 
arbitrary accuracy. This would replace the Brownian bridge correction usually 
combined with the equidistant scheme. We note that similar consideration 
motivated recently Gassiat et al. [81 who used Root's barrier hitting times for a 
1-dimensional Brownian motion in discretisation schemes. 


3.4 Implementation and numerical experiments 

For implementation purposes, note that the moving sphere scheme, or indeed 
the asymptotically efficient scheme in ll3.4b , have random time steps and we 
would typically overshoot the time horizon, i.e. have that > t. However 
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Figure 1: Reduction ratio rid) in J3.8t for the moving sphere scheme. 


in the case of moving sphere scheme this is easily corrected. Indeed, the bound 
A„,7 t" < means we can run the time steps until f - n™ > aG” and when this 
fails we may either chose to decrease G|5, or else continue with the equidistant 
scheme. This would involve at most 6 equidistant time steps in d = 1 and at 
most 4 for d>2. 

We report now a brief numerical study. We compare the moving sphere 
scheme with G = 1 and the standard equidistant scheme. Consider a two 
dimensional SDE 


/ tan((XHs) + X^(s))/2) ^ tan((XHs) - X^(s))/2) \ ^ 

1(1 + tan2((Xi(s) + X2(s))/2))2 (1 + tan2((Xi(s) - X2(s))/2))2 j 

+-^-7-dWi(s) +-^;-dW2(s) 

1 + tan2((Xi(s) + X2(s))/2) 1 + tan2((Xi(s) - X2(s))/2) 

/ tan((Xi(s) + X^(s))/2) tan((Xi(s) - X^(s))/2) \ ^ 

i(l + tan2((Xi(s) + X2(s))/2))2 (1 + tan2((Xi(s) - X2(s))/2))2 j 

+-^-r-dWi(s)-^;-dW2(s) 

1 + tan2((Xi(s) + X2(s))/2) 1 + tan2((Xi(s) - X^{s))/2) 


with X^(0) = X^(0) = 0. The solution admits an explicit expression 
X^(s) = arctan(W^(s)) + arctan(W^(s)), X^(s) = arctan(W^(s)) - arctan(W^(s)). 


We approximate the first four moments of the discretisation errors := 
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- X\l) and ;= X"'2(l) - X^{1) by the Monte Carlo with 1,000,000 
paths for both of the schemes. Normal random variables are generated by the 
Box-Miiller algorithm and an exponential variable is by - log(!J), where U is 
a uniform random variable. Uniform random variables are generated by the 
Merserme Twister algorithm. The Apple Mac mini computer with 2.6 GHz 
Intel Core i7 took 1 minute and 28 seconds for the equidistant scheme with 
n = 625. For the moving sphere scheme we took n = 435 which leads to an 
equivalent computational time of 1 minute and 26 second^ These two are 
therefore almost equivalent in terms of computation time. Table[T]and|2]report 
the first four moments of E"'^ and E"'^ respectively. 



E[E"A] 

E[|E”A|^] 

E[(E”A)J] 

E[|E"'T] 

equidistant (n = 625) 

8.1 xlQ-'’ 

0.00033 

1.1 xl0-» 

4.1 xlO-" 

moving sphere (n = 435) 

-4.7 xl0“‘’ 

0.00028 

1.9 xl0““ 

2.9 xl0“" 


Table 1: E”'^ = X"'i(l) - Xi(l) 



E[E""'] 

E[|E"'^|^] 

E[(E”'^)^] 

E[|E”'^|'‘] 

equidistant (n = 625) 

-1.1 xlO-^ 

0.00033 

-3.3 xl0-« 

4.1 xlO-" 

moving sphere (n = 435) 

-1.1 xlO-^ 

0.00028 

-7.8 xlO-^^ 

2.9 xlO-" 


Table 2: E”'^ = X"'2(l) - X2(l) 

From Table [T] and 121 we confirm that the moving sphere scheme provides a 
better accuracy without increasing computation time. 


4 Proofs 

4.1 Proof of Theorem 12.11 

Let em = inf{|x - yl;x e K,„, y e > 0, t„, = A t]™ and 

o’l„ = T]ni/\ infiw > 0;X„ i K,„ or Xjj ^ K„,+i}. 

Then 

P[tm < f] < P[ff", < f] < P[tm < f] + P[sup |!J„(S)| > Vwem]. 

0<s<f 

Since supg<j,<j |U"(s)| is tight by the assumption, for any e > 0, there exists X > 0 
such that 

limsupP[sup |U„(s)| > X] < e. 

n—>oo 0<s<f 

^ As described above, we follow the moving scheme algorithm until t- nil > and then we 
finish with 4 equidistant steps. This led to an average of 435.9 steps. 
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It follows that oJJj —> T„i, l{<j" <fi l{Tm<fi in probability and 


lim limsupP[a” < f] = 0. 

m^oo 


For any continuous bounded function cp on C([0, t], Rf), 


E[<p(U")] = +E[<p(U")l„»<„]. 


Therefore it suffices to show 


The coefficient / and its first derivatives are bounded and uniformly continuous 
on the compact sets K„+i. Therefore, for v < aJJ, A t, 


”''(U) r {/;'(X”(7t"(s))) - fl(X(s))j dW’is) 

j=0 *^0 

= ^L r { 4 (^"(^”(s))) - 4(^"(s)) + 4(^"(s)) - 4 (X(s))| dW^(s) 

/=0 *^0 

= -^E f '5fc/;'(^"(n"(s)))(X"''^(s)-X’''*^(7T”(s)))dW^(s) 

;,Jc 

+ f^k/;(X(s))(X'’'Hs)-X(s))dW^(s) + o,(l) 

j,k 

oo 

= - ^Y^Yj^kfi{X\nl))fl{X\ul)) I ' 

j,kj ^n=0 

+ £ f ^;c/;'(X(s))U''''^(s)dWi(s) + Op(l). 


(W'(s) - W‘(nl))dWi{s) 


Denote by V"'' the first of the two terms in the final expression above. Put 
X” = X”(n") and 

pn" Ac 

= (W'(s) - W'(7 t” _i))(W^(s) - W‘^(n:;_i))ds 

iJtI” ,Al 7 

m-1 

for 0 < b,c < d. Then V" = ..., P”'?) is a continuous semimartingale with 

quadratic covariation (P”'*, given by 

oo d d p 

-EEEE 

m=0 a=l b,c=0k,I=l 
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and 


oo p d 

dkfjiX^Jfl^iK) (W'(S) - W'«))ds. 

m=0 k=l 1=0 Jn^Av 

By Theorem IX.7.3 of Jacod and Shiryaev fTOl , if fhere exists a continuous 
process A = such that 

(V"i V'^'i) A'i W’) -» 0 (4.1) 


in probability as « —> oo for all i,], then V" converges !F-stably in law to a 
conditionally Gaussian martingale V = (V^,..VP) with {V’, Vi) = A'A 

We will argue below that (14.1b holds and that, using Z defined by (12.7b . the 
limit V is written as 


V{v) = 


d V pv 

= - E E d,fl{X{s))fl{X{s))AZ'’'\s). 

TTT Jo 


(4.2) 


a,b=l k=l 


The convergence of V" implies tightness of in C[0, f] by Theorem VI.4.18 of 

Jacod and Shiryaev IITOl . So any subsequence has a further subsequence which 
converges in law. Further it follows from (14.2b that the limit of the subsequence 
is uniquely determined by the SDE dl.Sb . Therefore U’\^„ itself must converges 
to Ll.At„ stably and we easily conclude. 

It remains to establish (14.lb . We do this in two steps. 

Step 1): We first show (W'G Wi) -> 0. 

By Ito's formula 



- W'«_i))ds = i(A” W')3 - p (W'(s) - W'(7T”_i))2dW'(s) 


for all 1 < Z < d. The conditional expectations of bofh terms in the right hand 
side are 0 by (12.5b . Further by (12.3b and (12.6b , 


K 

nY^E[\Ar„n’’\^\r:,_,] < 

m=l 

in probability and 


A 


N" 




m=l 


A 


K 


2^E[|A^7T«|2|!F;;_J 


m=l 


^ 0 
(4.3) 


E[|A" Wnr”_i] < CE[|A,„7 t”|3|^F;«_J, 

E[ (W'(s) - W‘{nl_,))^ds\r:_,] 

< lliminfE[|W'(M A nl) - W‘{nl_,)f\r"_,] < CE[|A„7T"nr,^_j] 
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for a constant C > 0. Then by Lemma A.2 of Fukasawa 13/ we obtain 


”■ (W'(S) - W'«,))ds ■ 

m=0 i:=l On” Ac 

in probability for 1 < Z < d. To treat the case I = 0, observe that 


0 




N," 


]E[|A„H”|2|r^_J < A £ ^ 0 

m=l ^ m=l 

by (12.4b and (12.6b . It follows then that (V"'', W^) —> 0 for all i,] again with the 
aid of Lemma A.2 of Fukasawa 0. 

Step 2): We show that (V”'', V"'^) converges and compute the limit 
By Ito's formula, (12.5b . we get 

for 1 < i>, c < d, where 6*’''^ is Kronecker's delta. The terms with = 0 or c = 0 
are negligible since 

K 

I ^ |A„ 


and 


,«|3 


0 


n; 

ny |A,„7t" 

m=l 

in probability, which follows from (12.6b and (14.3b by using Lemma A.2 of Fuka¬ 


sawa 0. Therefore, 


° a,i)=l A:,/=l *^0 


again by Lemma A.2 of Fukasawa |7|. This completes the proof. 

4.2 Proof of Lemma 

Proof: Let t be a stopping time with E[t] = a. Then W/^^ are uniformly 

integrable martingales and so, for any f > 0 by Jensen's inequality, 

d 

E[Q(T)rcAf] > Yj |E[W^'(T)|AcAt]l" = Q(t a f). 

/=! 

Therefore for (13.5b , it suffices to show 

E[Q(t)] > |^E[Tf 
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when T is a bounded stopping time. Let 


d 

S{v) = \Wi{vf. 

/=i 

Then 

d 

dQ(i7) = 4 2^(W^(s)fdW^'(s) + 6S(i;)di; 
y=i 

and 

d 

dS^(zy) = 2S(iy)dS(i7) + 4S(iy)di7 = 4S(iy)^ W’{v)dW’{v) + 2S{v){2 + d)dv. 

/=i 

It follows then 

E[Q(t)] = 6E[ S(i;)di;] = ^E[S(t)2] > ^E[S(T)f = 

The equality is attained if and only if S(t) is a constant, or equivalently, t is 
given by (13.6b . This completes the proof. 


References 

[1] S. Cambanis and Y.Hu (1996): Exact convergence rate of the Euler- 
Maruyama scheme, with application to sampling design. Stochastics Stochas¬ 
tics Rep. 59 211-240. 

[2] X. Chen, L. Cheng, J. Chadam and D. Sanders (2011): Existence and unique¬ 
ness of solutions to the inverse boundary crossing problem for diffusions, 
Ann. Appl. Probab. 21,1663-1693. 

[3] Z. Ciesielski and J. Taylor (1962): First passage times and sojourn times for 
Brownian motion in space and the exact Hausdorff measure of the sample 
path, Trans. Amer. Math. Soc. 103,434-450. 

[4] M. Davis, J. Obloj and P. Siorpaes (2015): Pathwise stochastic calculus and 
local times, in preparation. 

[5] M. Deaconu and S. Herrmarm (2013): Hitting time for Bessel processes - 
walk on moving spheres algorithm (WoMS), Ann. Appl. Probab. 23, 2259- 
2289. 

[6] H. Follmer (1981): Calcul d'lto sans probabilites. In Seminaire de Probabilites 
XV, volume 850 of Lecture Notes in Mathematics, pages 143-150. Springer. 

[7] M. Fukasawa (2011): Discretization error of stochastic integrals, Ann. Appl. 
Probab. 21,1436-1465. 


17 





[8] P. Gassiat, A. Mijatovic and H. Oberhauser (2015): An integral equation 
for Root's barrier and the generation of Brownian increments. Ann. Appl. 
Probah. 25,2039-2065. 

[9] J. Jacod and P. Protter (1998): Asymptotic error distributions for the Euler 
method for stochastic differential equations. Ann. Probah. 26 267-307. 

[10] J. Jacod and A. Shiryaev (2002): Limit Theorems for Stochastic Processes, 2nd 
ed. Springer-Verlag, Berlin. 

[11] PE. Kloeden and E. Platen (1992): Numerical solution of stochastic differential 
equations, Springer-Verlag, Berlin. 

[12] T.G. Kurtz and P. Protter (1991): Wong-Zakai corrections, random evolu¬ 
tions, and simulation schemes for SDEs. Stochastic Analysis 331-346. Aca¬ 
demic Press, Boston. 

[13] N. Landon (2013) : Stratgies de couverture presque optimale : 
thorie et applications. PhD thesis, Ecole Polytechnique, Paris, France, 
pastel.archives-ouvertes.fr/pastel-0IS788067. 

[14] G.N. Milstein, and M.V. Tretyakov (1999): Simulation of a space-time 
bounded diffusion, Ann. Appl. Probah. 9, 732-779. 

[15] T. Miiller-Granbach (2002): The optimal uniform approximation of sys¬ 
tems if stochastic differential equations, Ann. Appl. Probab. 12, 664-690. 

[16] N.J. Newton (1990). An efficient approximation for stochastic differential 
equations on the partition of symmetrical first passage times. Stochastics 
Stochastics Rep. 29 227-258. 

[17] N. Perkowski and D. Promel (2014): Local times for typical price paths 
and pathwise Tanaka formulas. Preprint, arXiv: 1405.4421. 

[18] V. Vovk (2012): Continuous-time trading and the emergence of probability. 
Fin. and Stoch., 16:561-609. 


18 


